Novel Mechanism for Discrete Scale Invariance in Sandpile Models 
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Numerical simulations and a mean-field analysis of a sandpile model of earthquake aftershocks in 
Id, 2d and 3d euclidean lattices determine that the average stress decays in a punctuated fashion 
after a main shock, with events occurring at characteristic times increasing as a geometrical series 
with a well-defined multiplicative factor which is a function of the stress corrosion exponent, the 
stress drop ratio and the degree of dissipation. These results are independent of the discrete nature 
of the lattice and stem from the interplay between the threshold dynamics and the power law stress 
relaxation. 



Discrete scale invariance (DSI) [ft] is the partial break- 
ing of continuous scale invariance [0 in which a system or 
an observable is invariant only under scaling ratios that 
are integer powers of a fundamental factor A. DSI leads 
to complex critical exponents (or dimensions), i.e. to log- 
periodic corrections to scaling, which reflect the existence 
of a discrete self-similar spectrum of characteristic scales 
decorating the usual scale-free power law behavior. 

Several mechanisms responsible for this partial break- 
ing of the continuous scale symmetry have been ex- 
pounded, which include build-in pre-existing hierarchy 
[3|, intermittent diffusion in discrete euclidean lattices 
\i\ and cascades of ultra-violet instabilities in growth 
processes and rupture ||. Other situations are less well 
understood but can be traced back to special technical 
properties such as the non-unitary structure of the un- 
derlying field theory describing the coarse-grained prop- 
erties of animals J(| and of quenched disordered spin sys- 
tems with long-range interactions [Q. Another example 
is the gravitational collapse leading to critical black hole 
formation described by a a system of partial differential 
equations possessing an asymptotic solution which can 
be understood from a renormalization group with a limit 
cycle having a single unstable mode 

Here, we present a novel scenario for DSI based on 
the interplay between the threshold dynamics charac- 
teristic of sandpile models and a scale-free relaxation 
process ||. Specifically, we study a conceptual sand- 
pile model (To) of earthquake aftershocks on a euclidean 
discrete ei-dimensional cubic lattice with L d cells and pe- 
riodic boundary conditions with d = 1,2,3. Each cell 
represents a region which is unloaded when an elemen- 
tary fault is activated. We neglect the tensorial nature of 
the stress field and consider an anti-plane driving config- 



uration in which loading and rupture are controlled by a 
single shear stress component l^(~a?). 

There are two distinct temporal phases. First, the 
stress is uniformly increased at a very slow rate on all 
cells to mimic the tectonic loading. Due to the rup- 
ture and loading rules described below, the system self- 
organises into a statistical stationary state, characterised 
by a power law distribution of event sizes JllJ|§- Once 
this statistical stationarity state is established, we freeze 
the loading and the aftershock sequence starts, mimick- 
ing the aftermath of a great earthquake. The second 
phase is characterised by the fact that the aftershocks 
are not driven by the tectonic loading but by relaxation 
processes as described below. 

An initial stress threshold B(~af) is assigned to each 
cell from a random uniform distribution in the interval 
Bq[1 — r, 1 + r]. We find similar results both for the 
annealed and quenched version of the model, in which 
either the thresholds are fixed or are resampled in the 
interval after each rupture. When the stress V(lt) in a 
cell at a? becomes larger or equal to B(!?), the stress is 
re-distributed according to the rules 
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v{i>)\t ei = y(^)|^ fOTC + n^)l bcforc Q-Jjpl . (2) 

Rule (||) applies to each of the 2d nearest neighbours 
(n.n.) of it carrying an initial stress ^("^)| b n foro which 
evolves into U("af)|^ cr . Because the toppling criterion 
(|]) depends only on V and not on its gradient, the order 
of site toppling commute jllj . 

7 is the relative stress drop, with 7=1 corresponding 
to a complete stress drop. (3 is the dissipation where 
1 — (3 quantifies the amount of stress drop transfered 
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to n. n. and is known as the seismic efficiency. For 
7 = 1, (0) and (|J) are identical to the rules used in the 
non-conservative sandpile model |l2] | , motivated from the 
coupling of blocks to a rigid upper driving plate in the 
Burridge-Knopoff model. Here, the dissipation accounts 
for the loss of stress and of stored elastic energy due to 
an earthquake under constant displacement conditions at 
the boundaries | Q . 

In the second relaxation phase, the loading stops and 
the thresholds decay in time according to the law 
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This model incorporates the mechanism of sub-critical 
crack growth and stress corrosion juj, which has been 
proposed as a possible delay mechanism for aftershocks 
|l5| , 16]|]] . In absence of loading, events are triggered each 
time the thresholds decay below the local stresses. When 
this occurs, the stress redistribution obeys ([!]) and 
The system eventually relaxes after a infinite time and in 
an intermittent manner to an equilibrium of zero stress on 
all elements of the lattice. A closed conservative spring- 
block system has also been found to relax in a infinite 
time and in an intermittent manner to the zero-stress 



equilibrium with self-organized critical behavior |18|. It 
is this complex relaxation that we study. It occurs via the 
triggering of what can be called aftershocks which exhibit 
remarkable properties. The results presented below are 
also found for a version of the model with continuous 
elasticity PJl9|| derived from Ref. |l7|] and are probably 
robust features of the general interplay between threshold 
dynamics and relaxation phenomena. 

We show here the simulations for 2d systems of 20 x 20 
elements (simulations have been performed with size up 
to 100 x 100 with no change of results) and in the an- 
nealed case where, after each toppling, thresholds are 
reassigned from the uniform distribution Bq[1 — r, 1 + r] 
with r — 0.75. The system is up-dated by finding the site 
closest to rupture and incrementing time, so that this site 
reaches its threshold. Once a site becomes unstable due 
to either loading (in the first phase) or by the decay of the 
threshold (in the second phase), stress is distributed to 
n.n. according to ([!]) and (||). The n.n. may also become 
unstable, releasing their stress and the process contin- 
ues until no further nodes are unstable, thus defining an 
event. During the rupture process, time is "freezed" to 
ensure a separation of time scales between rupture (fast) 
and loading/decay (slow). When no further sites are un- 
stable, loading or decay is continued until the next event, 
when time "freezes" again. Typically 3 • 10 3 to 15 • 10 3 
events were sampled in the decay regime. The results 
are robust with respect to heterogeneity level r, open or 
closed boundary conditions, and to the size and dimen- 
sion d = 1, 2, 3 of the lattice. 

Figure [l] shows the rate of aftershocks n (t) as a func- 
tion of time after the loading has ceased. The 1/t decay 



is in full agreement with Omori's law for real aftershock 
sequences |20) and is very robust over a wide range of pa- 
rameters, i.e. a > 0, 7 > and (3 < 1. Typical values for 
the earth are 10 < a < 100, 7 w 5 — 15% (stress drop) 
and (3 » 99% (corresponding to a seismic efficiency of 
1%). The distribution of event sizes also follows a power 
law (Gutenberg-Richter for aftershocks), but in contrast 
to Omori's law, the exponent continuously depends on 
the three parameters. Furthermore, for large dissipation 
/?, the power law extends only up to a maximum scale 
which decreases as (3 — > 1. 
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FIG. 1. Power law decay of the rate of aftershocks for a 
lattice of size 20 x 20 with a = 1, 7 = 1 and f3 = 0.5. The 
solid line is t~ . 



These results can be rationalized by the following mean 
field theory. From (||), we see that the time At needed 
for an isolated element to reach rupture is such that 
the threshold B (a? , to + At) decreases to the stress level 
V(lfr). The mean field argument simply assumes that we 
can extend this result over all elements of the lattice by 
replacing V( a?) by the average stress (V (t)) s accounting 
for the influence of the possible loading by n.n. This ap- 
proximation becomes better and better as the dissipation 
j3 increases and the dynamics of n.n. elements becomes 
increasingly uncoupled. This corresponds to a decreasing 
dependence on spatial inhomogeneities, which is exactly 
the underlying assumption of any mean-field approxima- 
tion. We get 
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where we have approximated Bq — (V)g by -Bo, since 
for large times the mean stress level becomes very low 
compared to the thresholds that are healed back to a 
typical value in the interval Bq[1 — r, 1 + r] after each 
event. Expression (Q) has the same form as obtained from 
a model of cracks undergoing sub-critical crack growth 
H . Over such a time interval, essentially one main event 
occurs on each site and, as a consequence, the average 
stress goes from (V) s to (1 — 7) (V)^, corresponding to 
a typical stress decrease 7 (V). We can thus write 
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whose solution is 



(^(i))^(So 2 /«7) 1/Q (t + c) 
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c is a constant determined from the initial value of the av- 
erage stress at the beginning of the aftershock relaxation 
sequence. 

To get Omori's law, we recognise that the rate n(t) of 
aftershocks is simply proportional to the rate with which 
the thresholds B (of, f) reach the stress level. n(t) is also 
proportional to 1/Ai. This yields n(t) oc dB (a?,t) /dt ~ 



with p = 1. Accord- 



(t+c)»> 

ing to this mean field theory, Omori's law is obtained 
with the universal exponent p = 1 independently of the 
value of the stress-corrosion exponent a, as long as it is 
positive. For a — 0, the average stress level decays ex- 
ponentially fast with time and the rate of aftershocks is 
constant. 

The mean field theory also provides a prediction of 
log-periodicity. The stress redistribution laws M) and 
(||) imply that, over a typical time At given by (Q), the 
average stress undergoes the change (V) s — > (V) s /fi, 
where 



= [f(d)(l - 7 ) + n cff7 (l - /?)] /f(d) . 
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/(ci) is a geometric factor counting the effective number 
of n.n. (f(d) — 2d in the large dissipation limit j3 — > 
1). The first contribution (1 — 7) in the r.h.s. of (J?]) 
is simply the initial stress minus the stress drop. The 
second contribution n c ff7(l — (3)/f(d) results from the 
number n e s ~ 1 of stress loading on a given element due 
to the earthquakes occurring on its neighbours. 

Each time the average stress is decreased by a factor 
fj,, we see from (^) that the time interval At is increased 
by a factor 



A = /i° 
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Since fi > 1, the total time is essentially dominated 
by the last time interval between the two last cycles. 
This allows us to write an approximate scaling rela- 
tion on the average stress (V) : (V (t)) s = [i (V (Xt)} s . 
Since the aftershock rate n(t) is proportional to (V)g 
(i = A, this leads to n(t) = X n(Xt), using (JsJ) . Look- 
ing for a power law solution n(t) ~ t~ p , we retrieve 
Omori's exponent p = 1. A more general solution is 
the power law t^ 1 multiplied by a periodic function of 
lni, n(t) — t~ 1 Pi(\nt/ In A), where P\ (x) is periodic with 
period one. Expanding this periodic function into its 
Fourier series gives n (t) = t^ 1 J2k=-oo a,kt l27lk / lnX , with 
a_fc = cifc, which defines the discrete spectrum of com- 
plex exponents pu — \ — ij^r- The leading correction to 
the power Omori's law gives the log-periodic expression 
n(t) = i(a + aicos(2^)). 




FIG. 2. Local exponent p(t) for a 
The window size is 100 events. 



To test this prediction, we find that the local ex- 
ponent p(t), defined by d\an{t) /dint = —p(t), gives 
the most sensitive measure of deviation from the 1/t 
Omori's law. We estimate p(t) by a maximum likeli- 
hood estimator in a running window ending at t ||. 
Defining the starting t and ending tjj times of a win- 
dow and the average (lnt^) of the logarithms of all 
N aftershock times within this window, the MLE is 
p(t) w 12(lnVHTy - (lnii))/(ln(t[//t)) 2 + 1), with a vari- 
ance a 2 w 12 (N - ly 1 {\ii{tu/t))- 2 . The estimation of 
p(t) is very robust over a large set of window sizes and 
have been tested thoroughly on synthetic Omori's laws 

Figure H shows the local exponent p (t) as a function of 
time t estimated using a window size of 100 events. Clear 
log-periodic oscillations around p ~ 1 can be identified 
with a (log-)frequency of « 0.12 giving a prefered scale 
factor A ps 3500 in reasonable agreement with the theo- 
retical value of 4000 calculated from . Increasing the 
window size with t or keeping a window with size fixed 
in time provides the same estimate. 
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FIG. 3. p- value as a function of time for the case of 
a — 25, 7 = 0.05, P = 0.99. The window size is 150 events. 



Figure (||) presents p(t) for the stress corrosion expo- 
nent a « 25, a value estimated from a set of adjacent 
time-delayed multiple events in western Japan |16|, for a 
stress drop 7 = 5% and a seismic efficiency of 1%, i.e. 
/3 = 0.99. The measured scaling factor is now A = 3.5 
while the mean field prediction is A = 3.6. 
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FIG. 4. Scaling factor A as a function of stress drop 7 for 
(3 = 0.999 and a = 1. 



The comparisons between the numerical simulations 
and the predictions of the mean field theory are presented 
in figures 4-6. The mean field theory, while not perfect, 
accounts well for most of the behavior. We also have 
checked that A depends on the space dimension d for 
d = 1, 2, 3 as predicted from 



log-periodicity does not rely on a pre-existing discrete 
structural hierarchy of faults but is dynamical and re- 
flects the existence of an approximately fixed stress drop 
together with the scale-free stress corrosion power law 
acting during inter-seismic phases. 

This study suggests to search for log-periodic signa- 
tures in real aftershock sequences, with the potential 
bonus that log-periodicity would constrain the stress 
drop ratio, an elusive quantity to estimate by direct seis- 
mic measurements. A systematic analysis on more than 
thirty large aftershock sequences found some indications 
H but was unable to conclude decisively yet J2TJ , due 
to a noise level which is comparable to the amplitude of 
the signal. Further studies are thus called for with better 
quality data. 

We are especially grateful to A. Johansen and L. 
Knopoff for very stimulating discussions. M.L. thanks 
the Physics department at UCLA for hospitality. 



FIG. 5. 
7 = 1. 
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Scaling factor A as a function of /3 for a = 1 and 
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FIG. 6. Scaling factor A as a function of a for j3 = 0.99 
and 7 = 1. 



Discrete scale invariance and its log-periodic signature 
is associated in this model to the threshold nature of the 
dynamics. The discrete scaling emerges from the fact 
that the thresholds are healed back to a value close to 
Bo after each rupture and then have to decay down to 
the current stress threshold to trigger the next rupture. 
When this occurs, the stress jumps to a smaller value by 
a finite amount and can also latter be reloaded again by 
active neighbours. It is fundamentally these finite jumps 
in the stress proportional to the current stress (which are 
thus scale-free) which are at the origin of discrete scale 
invariance and log-periodicity. This novel mechanism of 
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